package fem2.material.cohesive_law2;

import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;
import java.util.Map;
import java.util.TreeMap;

import math2.MathUtilities;
import math2.MutableDouble;
import fem2.Debugger;
import fem2.Element;
import fem2.Variable;
import fem2.element.crack_sda2.SDAElement;
import fem2.enu.VariableName;
import fem2.jlinalg.DiagonalMatrix;
import fem2.material.FluencyCriteria;

/**
 * Mode I cohesive law based on Sancho paper: An embedded crack model for finite
 * element analysis of concrete fracture.<br>
 * Note: this class works with constant opening only
 * 
 * @author hbui
 * 
 */
public class IsotropicDamageCohesiveLaw extends CohesiveLaw {

	private Map<Element, Boolean> isLoading = new TreeMap<Element, Boolean>();
	private FluencyCriteria softeningRule;

	public IsotropicDamageCohesiveLaw(double Gf, double ft, FluencyCriteria softeningRule) {
		super(Gf, ft);
		this.softeningRule = softeningRule;
	}

	/**
	 * compute value of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected double computeSofteningLaw(double x) {
		return softeningRule.computeFluencyValue(x);
	}

	/**
	 * compute value of derivative of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected double computeSofteningLawDerivative(double x) {
		return softeningRule.computeFluencyDerivative(x);
	}

	@Override
	public IVector getJumpTraction(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		IVector w = getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP);
		IVector t;

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);

		Debugger.watch("wtilde (t) = " + wtilde);

		if (wtilde < MathUtilities.tol6) {
			double ft = computeSofteningLaw(0);
			t = new ArrayVector(dim);
			t.set(0, ft);
			return t;
		}

		/*
		 * compute traction
		 */
		t = (IVector) w.clone();

		double alpha = getVariable(e).get(k).getScalarVariable(VariableName.ALPHA).getValue();
		double q = getVariable(e).get(k).getScalarVariable(VariableName.STRAIN_LIKE_VARIABLE)
				.getValue();

		Debugger.watch("alpha (t) = " + alpha);
		Debugger.watch("q (t) = " + q);

		BLAM.scale(q / alpha, t);

		return t;
	}

	@Override
	public IMatrix getJumpStiffness(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		IVector w = getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP);
		IMatrix C;

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);

		Debugger.watch("wtilde (C) = " + wtilde);

		if (wtilde < MathUtilities.tol6) {
			C = new DiagonalMatrix(dim);
			// BLAM.zero(C);
			return C;
		}

		double alpha = getVariable(e).get(k).getScalarVariable(VariableName.ALPHA).getValue();

		double alpha_old = getVariable(e).get(k).getScalarVariable(VariableName.ALPHA, -1);
		double q = getVariable(e).get(k).getScalarVariable(VariableName.STRAIN_LIKE_VARIABLE)
				.getValue();

		Debugger.watch("alpha (C) = " + alpha);
		Debugger.watch("q (C) = " + q);

		if (isLoading.get(e)) {

			Debugger.watch("loading");

			/*
			 * loading
			 */
			C = new Array2DMatrix(dim, dim);
			DiagonalMatrix.set(C, q / alpha);

			double H = computeSofteningLawDerivative(alpha_old);

			BLAM.multiply((H * alpha - q) / Math.pow(alpha, 3), w, w, C);

		} else {

			/*
			 * unloading
			 */
			C = new DiagonalMatrix(dim);
			DiagonalMatrix.set(C, q / alpha);

		}

		return C;
	}

	@Override
	public void Initialize(Element e) {
		int ngc = ((SDAElement) e).getNumGaussPointOnCrack();
		int dim = e.getMeshPart().getDimension();

		ArrayList<Variable> vars = new ArrayList<Variable>();

		for (int i = 0; i < ngc; i++) {
			Variable v = new Variable(1);

			v.registerVectorVariable(VariableName.DISCONTINUITY_JUMP, dim, false);

			v.registerScalarVariable(VariableName.STRAIN_LIKE_VARIABLE, true);
			v.getScalarVariable(VariableName.STRAIN_LIKE_VARIABLE).setValue(
					softeningRule.computeFluencyValue(0.0));

			v.registerScalarVariable(VariableName.ALPHA, true);

			v.save();

			vars.add(v);
		}

		setVariable(e, vars);
		isLoading.put(e, false);
	}

	@Override
	public void InitializeSolutionStep(Element e) {
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		int ngc = ((SDAElement) e).getNumGaussPointOnCrack();
		IVector w = ((SDAElement) e).getCrackOpeningParameter();
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);

		for (int k = 0; k < ngc; k++) {

			/*
			 * save the discontinuity jump
			 */
			BLAM.copy(w, getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP));

			/*
			 * check for loading/unloading
			 */
			MutableDouble q = getVariable(e).get(k).getScalarVariable(
					VariableName.STRAIN_LIKE_VARIABLE);
			double q_old = getVariable(e).get(k).getScalarVariable(
					VariableName.STRAIN_LIKE_VARIABLE, -1);
			MutableDouble alpha = getVariable(e).get(k).getScalarVariable(VariableName.ALPHA);
			double alpha_old = getVariable(e).get(k).getScalarVariable(VariableName.ALPHA, -1);

			if (wtilde > alpha_old) {
				isLoading.put(e, true);

				double dalpha = wtilde - alpha_old;
				alpha.setValue(alpha_old + dalpha);

				double H = computeSofteningLawDerivative(alpha_old);

				q.setValue(q_old + H * dalpha);

				Debugger.watch("H = " + H);
				Debugger.watch("wtilde = " + wtilde);
				Debugger.watch("alpha_old = " + alpha_old);
				Debugger.watch("dalpha = " + dalpha);
				Debugger.watch("q (old) = " + q_old);
				Debugger.watch("q (update) = " + q.getValue());
			} else {
				isLoading.put(e, false);

				alpha.setValue(alpha_old);
				q.setValue(q_old);
			}
		}

	}

	@Override
	public void FinalizeSolutionStep(Element e) {
		for (Variable v : getVariable(e)) {
			v.save();
		}
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

}
